Genetic structure shows the presence of small‐scale management units in a relict tree species

Abstract Identifying conservation units is crucial for the effective conservation of threatened species. Previous cases are almost exclusively based on large‐scale but coarse sampling for genetic structure analyses. Significant genetic structure can occur within a small range, and thus multiple conservation units may exist in narrowly distributed plants. However, small‐scale genetic structure is often overlooked in conservation planning especially for wind‐pollinated and wind‐dispersed trees, largely due to the absence of dense and elaborate sampling. In this study, we focused on a representative endangered relict plant, Metasequoia glyptostroboides. Using both nuclear microsatellites (nSSRs) and chloroplast DNA (cpDNA) fragments, we sampled across the narrow distribution range of this species and determined its conservation units by exploring its genetic structure and historical demography. cpDNA haplotypes were classified into two groups, but mixed in space, suggesting that the existent wild trees of M. glyptostroboides cannot be divided into different evolutionarily significant units. However, using nSSRs, we detected strong spatial genetic structure, with significant genetic differentiation and weak gene flow between the samples in the east of the species' distribution range and other samples. The divergence between the two nSSR groups was dated to the Last Glacial Maximum (c. 19.6 kya), suggesting that such spatial genetic structure has been maintained for a long term. Therefore, these two nSSR groups should be considered as different conservation units, that is, management units, to protect intergroup genetic variations, which is likely to be the outputs of local adaptation. Our findings highlight the necessity to reveal small‐scale genetic structure and population demography to improve the conservation strategies of evolutionary potential of endangered plants.


| INTRODUC TI ON
One of the most important mechanisms causing the endangerment of species is insufficient intraspecific genetic diversity to keep pace with the ongoing environmental changes (Fricke et al., 2022).Therefore, protecting genetic variations and the genetic distinctiveness is essential for the conservation of threatened species (Liu et al., 2020;Meirmans & Hedrick, 2011;Teixeira & Huber, 2021).To achieve this goal, a key step in conservation practices is to divide the extant populations into different conservation units, including evolutionarily significant units (ESUs), which are developed to provide a rational basis for recognizing evolutionary heritage to prioritize protected populations, and management units (MUs), which are genetically differentiated units that can be monitored independently (Funk et al., 2012;Moritz, 1994).
As significant genetic differentiation usually occurs over large regions, identification of conservation units was often performed at a broad scale (e.g., Olsen et al., 2014).However, populations can differentiate in a small range especially in mountainous regions due to local adaptation, founder effects by seeds stochastic colonization, and genetic drift in small and isolated populations (Frei et al., 2012;Steinfartz et al., 2007).Furthermore, given that distribution ranges of most threatened species are shrunken and that genetic differentiation within a small scale can be quickly formed in the presence of gene flow barriers (Vranckx et al., 2011), understanding the genetic structure to guide the design of smallscale conservation units of endangered species is thus becoming increasingly urgent.
When identifying conservation units, multiple types of molecular markers are required to identify ESUs and MUs (Funk et al., 2012).
To distinguish evolutionarily distinct taxa, ESUs should be reciprocal monophyly revealed by organelle DNA and show significant divergence of allele frequencies at nuclear loci (Moritz, 1994).For the convenience of conservation practices, groups of individuals with significant differentiations in allele frequency at nuclear or organelle DNA loci but not phylogenetic distinctiveness are also considered to be important, because they reflect intraspecific evolutionary potential and may represent functional differentiations, which is considered as MUs (Moritz, 1994;Olsen et al., 2014;Sandercock et al., 2022).Besides, another criterion of recognizing conservation units is that the observed strong genetic differentiations have been stable for the long term.However, this is seldomly tested in threatened species, mainly due to lack of knowledge about the demographic dynamics and the gene flow pattern.
Recently developed Bayesian approach has enabled the retrospection of historical demography using molecular markers and the inference of lineage split events and individual dispersal rates causing genetic differentiation (Beaumont, 2010;Setsuko et al., 2020), providing new insights into the delimitation of conservation units.
In contrast to the traditional estimation of gene flow indirectly according to genetic differentiations, which has low statistical power to delineate gene flow pattern across different historical periods (Palsbøll et al., 2007), this approach can assess the strength of both historical and recent gene flow among differentiated groups and subsequently explore the role of historical events such as the Last Glacial Maximum (LGM).Therefore, the origin of genetic differentiation and spatiotemporal pattern of gene flow are available, which are crucial to test whether genetic differentiation is stable in the extant populations.Nevertheless, information of both genetic structure and demographic dynamics is still absent for most threatened species, largely because genetic factors are less considered than demographic characteristics when designing conservation strategies (Liu et al., 2020).
Metasequoia glyptostroboides (Cupressaceae) is a well-known "living fossil" tree species as well as a flagship species.It was thought to have been extinct in the early Pleistocene but was rediscovered in the 1940s in Three Gorges Mountain Region (TGMR), China.Fossil records indicate that genus Metasequoia originated in the early Late Cretaceous (c. 100 million years ago) and once distributed over most areas of the Northern Hemisphere (LePage et al., 2005;Li et al., 2012;Yang, 1998).The distribution range of this genus began to retreat during the Eocene and Miocene (Yang, 1998).At present, wild M. glyptostroboides trees can only be found in a mountainous region of the border area of Hubei Province, Hunan Province, and Chongqing Municipality (c.800 km 2 ), China, with a census size of about 5500 (Figure 1a).Because of its narrow distribution range and small population size, M. glyptostroboides was listed as an endangered species of IUCN (Farjon, 2013), and the first-class of National Key Protected Wild Plants in China (http://www.forestry.gov.cn/).Our previous studies using random amplified polymorphic DNA markers detected a low genetic diversity in both its wild and planted populations (Chen et al., 2003;Li et al., 2005).Nevertheless, there are still no conservation measures concerning the protection of its genetic variations, largely because it is still unclear whether multiple conservation units exist within wild populations.
Complex landscape patterns often lead to population genetic structures through limiting gene flow and decreasing population size resulting in inbreeding and genetic drift (Love Stowell et al., 2017).
Given the high complexity of landscape features in its natural distribution range, we hypothesized that strong genetic structure and multiple management units presented in this species.To test this hypothesis, 467 trees were sampled across its natural distribution range and were analyzed using both nuclear microsatellites (nSSRs) and chloroplast DNA markers.We analyzed the genetic structure and history demography.Specifically, we answered the two questions: (1) are there genetically differentiated groups of sampled trees in space; and (2) if yes, when did they form and was there recent  Moritz (1994) and Palsbøll et al. (2007), to guide both in situ and ex situ conservation projects.

| Population information and tissue sampling
Metasequoia glyptostroboides is monoecious but mainly outcrossing.Both pollen and seed grains are dispersed by wind, with paternal inheritance of its chloroplast (Mogensen, 1996).Its typical habitats are humid forests and riparian zones.This species is distributed within a narrow and mountainous area with most wild trees (99% of wild trees) found in Lichuan County (Hubei Province, China), where the average annual rainfall is 1319 mm and the average temperature per month varies from 1.9 to 22.6°C (Li et al., 2021).In this area, M. glyptostroboides trees grow at elevations between 800 and 1400 m in two valleys (Xiaohe Valley in the west and Shiziba Valley in the east) separated by Hanpo mountain (the altitude of its peaks reaches 1400 m; Figure 1a; Tang et al., 2011).M. glyptostroboides is a dominant canopy species in local forests, and its adults can reach more than 30 m high with the diameters at breast height of adult trees of 30-350 cm (Li et al., 2021).In addition to Lichuan, scattered wild trees have been In this study, leaf tissues were collected from a total of 467 wild trees across its natural distribution range, including 443 trees from Lichuan, 21 trees from Shizhu, and three trees from Longshan (Figure 1a).The detailed sampling information is listed in Table S1.The collected fresh leaves were dried in silica gel for DNA extraction.

| Genetic diversity and phylogenetic history based on cpDNA data
To estimate cpDNA diversity, we calculated number of haplotypes (H), haplotype diversity (h), and nucleotide diversity (π) using DNASP v. 6.12.01 (Rozas et al., 2017).Deviations from neutrality of sequence variations were tested by Tajima's D statistic and by Fu and Li's D* and F* statistics using DNASP v6.12.01.
To test whether there were multiple genetically differentiated groups within the sampled trees, a haplotype network was constructed using statistical parsimony as implemented in TCS v.1.2(Clement et al., 2000).Ambiguities in the network were resolved using the criteria and procedures delineated by Crandall and Templeton (1993).
Furthermore, we constructed the phylogeny and estimated divergence time between clades constructing Bayesian phylogenetic trees using BEAST v.1.7.5 (Drummond et al., 2012).The sequences of Cryptomeria japonica (Taxodiaceae; GenBank accession number: AP009377), Taiwania cryptomerioides (Taxodiaceae; GenBank accession number: NC_016065), and Cunninghamia lanceolata (Taxodiaceae; GenBank accession number: NC_021437.1) were used as outgroups.The best model for nucleotide substitutions (the HKY + I + G model) was selected by the Akaike information criterion using jModeltest 2.1.4(Darriba et al., 2012).Given that the phylogenetic tree was constructed for gymnosperms, a plant group with high probability of extinct species, a birth-death speciation model with incomplete sampling was used as tree prior with estimated base frequencies (Stadler, 2010).Calibration points estimated from a molecular clock (Mao et al., 2012) were treated as priors on clade ages with a normal distribution.We used three molecular clock models (strict clock model, uncorrelated relaxed lognormal clock model, and uncorrelated relaxed exponential clock model).Three independent analyses were run for 30 million generations, sampling trees at every 3000 generations.Convergence to the posterior distributions of divergence times and parameter estimates was examined in TRACER v.1.5(Rambaut & Drummond, 2007), and the burn-in was set to be 10%.Based on Bayes factors calculated by TRACER, uncorrelated relaxed exponential clock was selected as the optimal clock model.Using this model, another two independent analyses were run.A maximumcredibility tree was calculated using TREEANNOTATOR v.1.7.4 (Drummond et al., 2012) after removing trees sampled during the burn-in phase.
To estimate the divergence time, we chose three molecular clock calibration points from Mao et al. (2012)  Genetic variation structure was estimated by analysis of molecular variance (AMOVA, Excoffier et al., 1992) on clade levels in GENALEX version 6.4 (Peakall & Smouse, 2006).Pairwise ′ ST was then estimated by φ ST /φ ST(max) (Meirmans, 2006) between any two clades to infer the strengthen of pollen-mediated gene flow, using the AMOVA approach in GENALEX.

| Genetic structure and genetic diversity based on nuclear microsatellite data
Genetic diversity of eight SSR loci was estimated in GENALEX software using the following statistics: number of alleles (N a ); number of effective alleles (N e ); observed heterozygosity (H O ); expected heterozygosity (H E ; Nei, 1978); and inbreeding coefficient (F IS ).We allocated all sampled individuals into 18 populations according to the locations of sampling sites to perform the linkage disequilibrium test and null allele frequency analysis (Table S2).Departure from Hardy-Weinberg equilibrium and linkage disequilibrium were examined using the exact tests in GENEPOP v.4.2.2 (Rousset, 2008) with sequential Bonferroni adjustment (Rice, 1989).The presence of null alleles was estimated using the software Micro-Checker v.2.2. 3 (van Oosterhout et al., 2004).
Bayesian clustering was applied to characterize population structure and estimate genetic admixture in all individuals using the software STRUCTURE version 2.3 (Pritchard et al., 2000).We ran 20 replicates of 10,000 Markov Chain Monte Carlo (MCMC) steps after 10,000 burn-in steps for each of 1-15 genetic clusters (K).Admixture model and correlated allelic frequencies were used to analyze the data set.We used the ΔK, which measures the rate of change in the log probability of the data between successive K values (Evanno et al., 2005), to determine the optimal value of K.For the selected K value, clustering analysis was performed setting 1000,000 MCMC steps following 100,000 burn-in steps.The output was visualized using the program DISTRUCT (Rosenberg, 2004).
Genetic diversity and differentiation between each pair of groups were estimated by Wright's F ST using GENALEX, as it is an appropriate index when using molecular markers with high mutation rates such as microsatellites (Whitlock, 2011).Moreover, we also calculated F ′ ST (i.e., F ST /F ST(max) ) to obtain standardized measures of genetic differentiation using AMOVA (Meirmans, 2006) in GENALEX.AMOVA was also performed to quantify genetic variations among groups, as implemented in GENALEX program.
Principal coordinate analysis (PCoA) was used to visualize genetic divergence among samples in a multidimensional space based on distance matrix with data standardization, using GENALEX program.
To test the phylogeographic structure, we compared the global empirical R ST and permuted R ST (pR ST ) using 10,000 permutations across allele size using the program SPAGeDi (Hardy & Vekemans, 2002).A significantly larger R ST than pR ST suggests that phylogenetically close alleles are more aggregated in a population than expected by random situations (Hardy et al., 2003;Hmeljevski et al., 2017).
Recent migration within the last several generations among each pair of groups was estimated using BayesAss v.3.0.4 (Wilson & Rannala, 2003).Markov chain Monte Carlo analysis was done with 1000,000 iterations, 100,000 burn-in length, and 1000 samples.To examine the convergence, we used TRACER v.1.5(Rambaut & Drummond, 2007) to analyze the trace file for the log probability of runs.We also estimate historical migration rates using Bayesian inference based on structured coalescent by Migrate-n version 4.4.3(Beerli et al., 2019).A Brownian motion model of mutation was used under Bayesian inference.We used three long chains 5,000,000 iterations in length and four static heating chains at temperatures 1.00, 1.50, 3.00, and 1000,000.

| Demographic history estimation
The demographic history was also detected with an approximate Bayesian computation (ABC) procedure (Beaumont, 2010) following the Setsuko et al. (2020) and Jiang et al. (2021) approach.We estimated the population dynamics, population divergence, and historical migration based on all sampled trees using nSSR loci.
In the analysis of population size change, we calculated the average and standard deviation of the number of alleles, heterozygosity, and allele size range using ARLSUMSTAT v.3.6 (Excoffier & Lischer, 2010).For each genetic group, we compared the follow-  S3.
The simulations of three models for nSSR data were performed by fastsimcoal2 v.2.6.0.3 (Excoffier & Foll, 2011).We used the generalized stepwise mutation model (GSM), setting the mutation rate per generation (μ) as 5 × 10 −4 with the prior distribution for each locus drawn randomly from gamma distribution with shape and rate parameters (Setsuko et al., 2020) and the geometric parameter (P GSM ) ranging from 0 to 1 for eight independent loci.All priors were generated by R v.4.0.3 (R Core Team, 2018).The three population size change models were simulated 100,000 datasets and summary statistics were calculated using ARLSUMSTAT.We then compared these three models by the ABC random forest approach using the R package ABCRF v.1.7 (Pudlo et al., 2016), with the number of trees in random forest setting as 1000.The classification error and posterior probability for each model were calculated to determine the best model.The goodness of fit for each scenario and priors was measured by the discrepancy between 1000 simulated pseudo-observed datasets and our experimental data.
Then three different population divergence models were compared (Figure 2).In population divergence analysis, the average total allele size range and G ST(b) between the two clusters were calculated using ARLSUMSTAT, setting the other parameters using N CUR1 and N CUR2 were fixed according to the results in single population size model (Table S3).In all models, the period of migration was set from T to the present.Prior distribution of divergence time (T), ancestral effective population size (N ANC ), and migration rate (m) were drawn from uniform distribution from 1 to 10 4 , from 10 2 to 10 4 , and from 10 −5 to 0.3, respectively.Simulations and model choice were conducted using the same approaches in the analysis of population dynamics.For the optimal population divergence model, we conducted 10 5 simulations, with posterior distribution of parameters estimated using neural network regression by the R package ABC v.2.1 (Csilléry et al., 2012).The posterior mode and 95% highest posterior density were assessed using the density function and the R package coda (Plummer et al., 2006).The goodness of fit for the optimal model was evaluated using the difference between 1000 simulated pseudo-observed datasets and our experimental data.

| Delimitation of conservation units
Based on our results, we tried to delimit different types of conservation units for M. glyptostroboides using the following standards.
Evolutionarily significant units (ESUs) were identified according to the definition from Moritz (1994), that is, within-species lineages that were monophyletic for plastid DNA alleles with significant divergence of allele frequencies at nDNA loci.Management units (MUs) are different groups of populations with significant divergence of allele frequencies at nDNA or plastid DNA and of <10% of the dispersal rate of individuals between groups (Moritz, 1994;Palsbøll et al., 2007).
After sequence alignment and removal of ambiguous regions, we obtained 2827 aligned nucleotides for all M. glyptostroboides haplotypes and three taxa (C.japonica, T. cryptomerioides, and C. lanceolata) for constructing phylogenetic trees.Using the Bayesian optimization approach, we found that all haplotypes of M. glyptostroboides were divided into two clades with a high posterior probability (1.00) (Figure 3), and the composition of haplotypes in each clade (centered on haplotypes 1 and 7, and centered on haplotypes 9 and 10, respectively) exactly matched with each group detected in the haplotype network (Figure 3).In addition, the phylogenetic relationships among haplotypes within each clade were poorly resolved (posterior probabilities <0.5 for all within-clade nodes).These results showed the presence of two genetically differentiated clades in M. glyptostroboides based on cpDNA data.The divergence time between the two clades was estimated to be c.9.4 Mya (95% HPD: 1.4-20.5Mya) in Neogene (Figure 3).
Moderate genetic differentiation was found between the clades of sampled trees ( ′ ST = 0.241), and AMOVA did not revealed a high intergroup differentiation (F CT = 0.022, p = .003,Table S5).However, the sampled individuals assigned to the two genetic clades scattered in the study region and did not show a genetic structure in space (Figure S1).

| Genetic diversity and genetic structure based on nuclear microsatellite data
Using the eight SSRs for all the 467 sampled trees, we detected 50 alleles in total, and H O and H E were 0.44 and 0.44, respectively (Table 2).The inbreeding coefficient per locus ranged from −0.648 to 0.541 with significantly high or low values in six loci, but the overall F IS across all loci was only 0.003 (SE ± 0.041), showing a low level of inbreeding (Table 2).We did not detect significant linkage disequilibrium in any pair of loci.In addition, there was no solid evidence for the presence of null alleles (Table S6).We found a significant departure from Hardy-Weinberg equilibrium in loci Mg23, Mg61, and Mg76, due to a heterozygote deficit reflected by the positive F IS values (Table 2).The genetic clustering analysis using STRUCTURE showed that when K = 2, most sampled trees located in different valleys were assigned into distinct genetic clusters (Figure 1c), forming the east and the west groups mainly segregated by the mountain lying between the two valleys (Figure 1a).Further division of sampled trees only occurred in the west group at K = 3 and 4 (Figure 1c).
PCoA also clustered all sampled trees into an east and a west group, with the same assignment results as those in STRUCTURE analysis at K = 2 (Figure 1b).Therefore, these results supported the presence of two groups, despite that the optimal K was 3 according to the results of ΔK analysis (Figure S2).Note that the assignment of sampled trees based on nSSRs was inconsistent with the two cpDNA-based genetic clades, which did not display apparent genetic structure in space.
We found high genetic differentiation between the east and the west groups (F ST (0.112) and F ′ ST (0.213)), and AMOVA also showed a significant genetic differentiation between the two groups explaining 18.34% of total genetic variance (Table S5).A phylogeographic structure was detected by a significant difference between the global estimates of R ST and permuted R ST (R ST = 0.2268, pR ST = 0.0006, p < .001).
Moreover, we detected weak but significant recent gene flow between the two groups, as of the lower boundaries of 95% confidence interval of migrate rates per generation (m) revealed by BayesAss were larger than 0, with the m from the west group to the east one being 0.0077 (95% CI: 0.0010-0.0208)and the m in the reverse direction being 0.0085 (95% CI: 0.0012-0.0192).A similar pattern of historical gene flow was also found using Migrate-n, showing very low values of m (from the west to the east group: 0.0002 (95% CI: 0.0001-0.0003)and the reverse direction: 0.0004 (95% CI: 0.0002-0.0005)).

| Demographic history and conservation units
Among the three population size change models, the standard neutral model (SNM) was the optimal model, which had the highest posterior probabilities for both groups.Average classification error rates of compared models were 0.295 and 0.315 of the west and east groups, respectively (Table S7).This result suggested that the population sizes of both groups were stable (Figure 2).In population divergence analysis, the optimal model was the divergence model with gene flow from the east group to the west one (DVM1) with a posterior probability of 0.502 (Table S7).In addition, for the optimal model in each analysis,  TA B L E 2 Genetic diversity calculated for eight loci and two groups from nuclear microsatellite data of 467 Metasequoia glyptostroboides sampled trees.
1975 for the west group (N CUR1 ), the east group (N CUR2 ), and their ancestral population (N ANC ), respectively (Table 3).The posterior mode of divergence time T was assessed to be 783 generations ago (19,575 years ago with the generation length of 25 years; Table 3).
Moreover, the posterior mode of the number of migrants per generation (Nm 21 ) from the east group to the west one was very low (1.738), with a migration rate of 0.0013 (Table 3).
Because there was no clear structure of the two cpDNA-based clades, the wild M. glyptostroboides trees could only be considered as one ESU.However, using nSSRs, based on the high genetic differentiation of allele frequencies showed by STRUCTURE and PCoA (Figure 1), and weak gene flow detected by BayesAss and demography study between the two groups for nDNA loci, the wild M. glyptostroboides trees located in the east and west of the study area should be considered as two MUs.

| DISCUSS ION
Identifying conservation units is of great importance for the conservation of threatened species (Moritz, 1994).We identified two MUs in wild M. glyptostroboides trees using nSSR data, which locate in the east and west of our study area, with significant genetic differentiation and weak gene flow.However, it cannot be divided into different ESUs based on cpDNA information.Given that genetic factors have not yet been considered in the conservation of this species (Li et al., 2005), the determination of conservation units will greatly improve the conservation strategy of M. glyptostroboides.
The phylogenetic tree showed two clades in M. glyptostroboides, but the haplotypes mixed in space, supporting the presence of only one ESU.M. glyptostroboides was widely distributed across the northern hemisphere in the Paleocene (Lepage et al., 2005), and began to contract in the late Miocene and Pliocene probably due to the cool temperate and arid climate developed (Lepage et al., 2005;Zhang & Wang, 2015).Fossil records indicate that Metasequoia distribution continued to contract, and the remnant populations in China and Japan were under the control of warm temperate and moisture climate, until disappeared from Japan and other locations in China except the TGMR during the early Pleistocene (Liu et al., 2007).
Thus, those haplotypes from different regions were admixed during the contractions and survived in the sole refugium.
Spatial Bayesian methods can overestimate the number of discrete clusters when faced with isolation by distance (IBD) pattern (Frantz et al., 2009).We thought two discrete groups could be identified across M. glyptostroboides distribution based on nSSR data by STRUCTURE and PCA, despite three clusters indicated by STRUCTURE without observed IBD and obvious biological explanation (Frantz et al., 2009).The refugium population of M. glyptostroboides was estimated to be formed at the early Pleistocene (Liu et al., 2007).This is consistent with the factor that the TGMR was one of the most important refugia for several ancient "tertiary relics" throughout the Quaternary, such as Ginkgo biloba (Hu, 1980;Shen et al., 2005), probably because effects of glaciations were far weaker than those in Europe and North America (Hu, 1980).Thus, we can confirm that the genetic structure (happened during the LGM) was formed in the current refugium.Moreover, the formation of two groups may be associated with the mountain separating two valleys.During the last glaciation, the mountain ridge might not have been an appropriate habitat (probably due to cold and dry climate) for M. glyptostroboides, which might have retreated to the areas of lower elevations, creating a long-term discontinuous distribution.Though M. glyptostroboides is a wind-pollinated and -dispersed species, it is still very difficult to disperse its propagules upwards over the ridge without any stepping stone conspecific trees or habitats.This is supported by the low but significant recent and ancient migration rates between the two groups.Thus, the mountain ridge was highly likely to act as the gene flow barrier, leading to the significant, though small-scale, differentiation in nuclear DNA.Nevertheless, the warming climate after LGM may assist the return of M. glyptostroboides to the high elevations of the mountains and weaken the effects of gene flow barrier.This is probably why a few sampled trees assigned into the dominant genetic cluster in the western valley appeared in the eastern valley.
Difference spatial genetic patterns were revealed using two types of genetic markers.This is probably because chloroplast DNA and nSSR markers are designed to explore the genetic variations generated at distinct time scales.Chloroplast DNA is paternalinherited in gymnosperms, of which genetic variations mainly result from accumulated mutations with little or no crossing over (Mc-Cauley, 1995), and thus are able to probe the impacts of historical events.However, nSSR markers are expected to detect the recently formed genetic structure, due to the far higher mutation rate in nSSR regions than chloroplast DNA and frequent recombination of nSSR loci (Petit et al., 2005).Species with limited distribution ranges and small population sizes often attract extensive attention of conservation because they face high risk of extinction (Pimm et al., 2014).To preserve the adaptative potential of an ESU in the face of rapidly changing environmental pressures, the total genetic variation harbored by all wild M. glyptostroboides trees should be maintained integrally (Li et al., 2020).Nevertheless, conservation practices should be conducted independently in each MU to protect the evolutionary potential of this relict species (Teixeira & Huber, 2021), as well as to achieve pathogen and pest prevention (Olsen et al., 2014).Thus, human-mediated migration between MUs should be avoided, leading to the question that whether the two isolated MUs can persist.Generally, an effective population size (EPS) of 500 is considered as a threshold for long-term persistence of a threatened species (Liu et al., 2020).Given the EPS of either MUs are larger than 500, such isolation is not likely to detrimentally impact MUs in the genetic aspect.In each MU, identification of wild individuals (e.g., trees containing rare alleles) for conservation priority must be first carried out to determine the core germplasm resources of M. glyptostroboides (Lu et al., 2007).Then, we can reinforce gene flow across different areas within each MU by facilitating developments of seeds and seedlings (Li et al., 2021), to maintain high levels of genetic diversity and low risk of inbreeding depression.
Moreover, human-assisted regeneration strategy like weeding is also necessary to improve the survival and growth of M. glyptostroboides seedlings and saplings (Li et al., 2021), for ensuring the effects of reinforced gene flow and enhancing the effective population size.
As to the ex situ conservation of M. glyptostroboides, although numerous planted populations have been established around the world, all of them are suffering from extremely low genetic diversity (Li et al., 2005), threatening their persistence in the presence of rapid environmental changes and pests/pathogens (Li et al., 2012).Therefore, it is essential to elevate the genetic variation in planted populations by collecting seeds from multiple trees in the core germplasm resources.However, it is worth noting that germplasms from different MUs should be introduced into distinct planted populations rather than introducing mixed germplasms (Li et al., 2005).Moreover, migration from planted population to the wild population must be prevented though there were many plantations nearby the distribution area of wild trees, in case of the Nowadays, drastic decline in population sizes and deprival of habitats by global changes and human disturbances have been reported in an increasing number of species (Díaz et al., 2019;Pimm et al., 2014).However, insights into small-scale genetic structure and demographic dynamics are almost overlooked for threatened species with small distribution ranges, leading to less effective conservation measures.Our study provides an example that conservation units divided by the locally originated genetic structure can even exist within a very small range, highlighting the necessity to conduct small-scale research when evaluating the existing genetic diversity for threatened species.
differentiation.We will delimit ESUs (by chloroplast DNA sequences) and MUs (by nSSRs) according to the definition proposed by

F
Sampling locations and genetic clustering of wild Metasequoia glyptostroboides trees using nSSR data.(a) The upper left figure shows the location of distributed area.The lower left figure shows the sampling region in the study.Blue backgrounds show the distribution range of the species.Red frame represents the region showed in the right figure, and dots at the north and south represent trees in Moudao town, Lichuan County, Hubei Province, and Luota town, Longshan County, Hunan Province, respectively.The right figure shows the main sampling area.Dots represent the locations of trees with different clusters by STRUCTURE (K = 2) using nSSR data, showing the topography and river network of the main distribution area.Blue lines represent the rivers along the mountainous valleys.Red dashed line indicates the division of the east and west groups according to the result of STRUCTURE.(b) Principal coordinates analysis (PCoA) based on nSSR genetic distances of individuals.The two genetic clusters inferred by STRUCTURE are represented by different colors.(c) Each individual is shown as a vertical line divided into segments representing the estimated membership proportions in the two, three, and four ancestral genetic clusters using STRUCTURE of nSSR data.Individuals within each cluster are arranged according to estimated cluster membership proportions.found in Shizhu (Chongqing Municipality) and Longshan (Hunan Province).
ing three models (1) the standard neutral model (SNM) assuming a constant population size (N) over time; (2) the population growth model (PGM) where the ancestral effective population size N a is assumed to increase exponentially (experiencing t generations ago) with a growth rate of G, that is, N a = N × exp (G × t), to reach its current effective size (N); and (3) the size reduction model (SRM) which describes a population with an ancestral effective population size N b experienced population decline from instantaneously t 1 generations ago to reach its current effective size (N, N < N b ).Details for prior distribution of parameters for each model are shown in Table the same statistics of the best population dynamics model.The divergence model assumed the divergence between different genetic groups from their ancestral population with effective population size N ANC at time T.Both unidirectional and bidirectional migration rates were set between any two genetic groups, forming three population divergence models.Each model has basic structural parameters: effective population size of a genetic group (such as N CUR1 and N CUR2 ), ancestral effective population size (N ANC , N CUR < N ANC ), divergence time (T), and migration rate (m).

F
Demographic models compared in approximate Bayesian computation.(a) SNM (standard neutral model): the population size was constant at all the time.N 1 and N 2 were the current effective population sizes of the west and east groups, separately.PGM (population growth model): the population expanded at the time T. N a1 and N a2 were the effective population sizes of the west and east ancestral populations, respectively.SRM (size reduction model): the population size contracted at the time t 1 .N b1 and N b2 were the effective population sizes of the west and east ancestral populations, respectively.t 1 = divergence time-scaled by generation time; N # = effective population size for two groups.(b) Divergence model (DVM) and dispersal model (DPM) at T generations ago with different gene flow.Arrows indicate the presence of gene flow.| 7 of 13 LI et al.
most observed summary statistics were close to the simulated values based on parameter values generated from the posterior distributions (Figures S3 and S4), suggesting a high goodness of fitting.Under these two selected models, the posterior mode values of effective population sizes were estimated to be 1327, 588, and F I G U R E 3 Chronogram for 47 haplotypes of M. glyptostroboides and three outgroups (Cryptomeria japonica, Taiwania cryptomerioides, and Cunninghamia lanceolata) based on chloroplast DNA fragments.Upper numbers represent divergent times (Mya, millions of years ago) with the 95% confidence interval of the mean time to the common ancestor in brackets, and lower numbers represent posterior probabilities.Phylogeny and divergent time are estimated through Bayesian methods in BEAST.Tip colors indicate the haplogroups as inferred by the phylogeny.Inserted small picture at the center indicates network of cpDNA haplotypes of M. glyptostroboides.The haplotype colors refer to the phylogeny analyses.The area of each circle is proportional to the frequency of each haplotype.Each solid line represents one mutational step that interconnects two haplotypes.The small hollow circles indicate ghost haplotypes.
loss of genetic distinctiveness due to genetic admixture.Although our findings can substantially improve current conservation strategy, future studies are expected to test whether differentiations in adaptation to environment conditions exist between MUs for the determination of suitable habitats for ex situ conservation and the prediction of the fate of M. glyptostroboides under the ongoing global changes.
IS , inbreeding coefficient; H E , expected heterozygosity; H O , observed heterozygosity; N a , number of alleles; N e , number of effective alleles.
Parameter estimates for the best demographic models based on approximate Bayesian computation of the west, east, and both groups.GSM , mean value of parameter in generalized stepwise nutation model among eight loci; N ANC , the ancestral effective population size; N CUR1 , current effective population size of west group; N CUR2 , effective population size of east group; Nm 21 , number of migrants per generation from east group to west one; Shape, a parameter of geometric distribution related to variance among eight loci; T, divergence time of west and east groups from their ancestral population.
TA B L E 3